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ABSTRACT 

The velocity dispersion of nearby stars in the Galactic disc are well known to 
increase substantially with age; this is the so-called Age- Velocity relation, and is in- 
terpreted as a "heating" of the disc as a function of time. We have studied the heating 
of the Galactic stellar disc due to giant molecular clouds and halo black holes, via 
simulations of the orbits of tracer stars embedded in a patch of the local Galactic disc. 
We examine a range of masses and number densities of the giant molecular cloud and 
halo black hole perturbers. The heating of the stellar disc in the simulations is fit with 
a simple power law of the tr oc where a is the velocity dispersion of the tracer stars 
as a function of time, t. We also fit this form to the best determinations of the in- 
crease in the velocity dispersion as a function of time as derived from stars in the solar 
neighbourhood for which ages can be reliably assigned. Observationally, a is found to 
lie in the range 0.3 to 0.6, i.e. it remains poorly constrained and its determination is 
probably still dominated by systematic errors. Better constrained observationally is 
the ratio of the velocity dispersions of the stars in the vertical z and horizontal x (i.e. 
toward the Galactic center) directions, being crz/<Jx = 0.5 ± 0.1. 

For the heating of the stellar disc due to giant molecular clouds (GMCs) we derive 
a heating a oc t'^''^^, which differs somewhat from early (analytic) studies in which 
a cx t'^-^^. This confirms the well known results that there are insufficient GMCs to 
heat the Galactic disc appropriately. A range of dark halo black hole scenarios are 
verified to heat the stellar disc as cr (x t'^^^ (as expected from analytical studies), and 
give az/(Jx in the range 0.5 to 0.6, which is consistent with observations. Black holes 
with a mass of 10'' Mq are our favoured disc heaters, although they are only marginally 
consistent with observations. Simulations featuring a combination of giant molecular 
clouds and halo black holes can explain the observed heating of the stellar disc, but 
since other perturbing mechanisms, such as spiral arms, are yet to be included, we 
regard this solution as ad hoc. 

Key words: Galaxy: kinematics and dynamics - (Galaxy:) solar neighbourhood - 
Galaxy: velocity dispersion - N-body simulations 



1 INTRODUCTION 



The rai dom motions of stars near the Sun are well known to 



mcreas* with stellar age, and this ettect is known as the disc 
Age Velocity Halation (hereatter AVKj. Broadly speaking, 
the total velocity dispersion of stars rises from ~20 kms~^ 
at the lowest measurable ages to « 60-80 kms~^ at an age 
of ~ 10-12 Gyr. The "heating" of the orbits is initially rapid 
but levels off after some Gyr. The heating is not uniform in 
each of the three velocity components; the velocity disper- 
sion rises to a higher value in the Galactic plane than in the 
vertical direction. 

A number of mechanisms have been proposed to explain 
the AVR, in particular the gravitational perturbative effect 



on stellar orbits by giant molecular clouds (GMCs) in the 
Galaxy's gas layer ( Spitzer & Schwarzschild 1951; Spitzer & 



Schwarzschild 1953), by spiral arms in the disc (Barbanis & 
Woltjer 1967[), or by massive black ho les in the Galactic dark 
halo (BHs) (Lacey & Ostriker 1985). Dramatic disc heat- 



ing als o takes place when satellite galaxies fall onto galactic 
discs (Quinn fc al. 1993), and such events may well leave 



an abrupt feature in the AVR, but if the rain of such satel- 
l ites is relatively gentle th ey might also heat discs smoothly 
( [Velazquez fc White 199^ . 

Spitzer and Schwarzschild (1 



1951 



1953 ) first argued that 



the observed growth of the stellar velocity dispersion could 
be explained by encounters of disc stars with GMCs with 
masses of Mgmc ~ IO'^Mq. The main difficulty with this 
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scenario is that the observed number of GMCs modify stel- 
lar orbits in a manner inconsistent with observations. Firstly, 
they are too few to reproduce the observed heating, and sec- 
ondly, too much heating is produced in the vertical direction 
relative to the disc plane ( |Lacey 1984 ). However, GMCs do 
exist so they certainly provide a mechanism for redirecting 
orbits out of the Galactic plane, whatever is responsible for 
the heating. Heating due to transient spiral arms suffers the 
opposite problem that the amo unt of vertical heating is too 
low compared to observations ( Carlberg 1987 ). 

This naturally leads to the idea that the observed heat- 
ing could be explained by the combined effect of the tran- 
sient spiral arms and giant molecular clouds. In this model 
spiral density waves would act as main source of heating 
and GMCs would delect the stellar orbits out of the Galac- 
tic plane, thus creati ng the required velocity dispersion in 
the vertical direction (Carlberg 1987). However, these calcu- 
lations are rather complicated: e.g. the relative effectiveness 
of spiral and GMC heating has to modelled with an em- 
pirical parameter. The rate at which spiral features heat a 
stellar disc depends on the potential's spatial pattern and 
on the time variability of the pattern. In principle, spatial 
stuctures of galaxies can be determined from photometry of 
real galaxies and temporal information about spiral stucture 
can be obtained from numerical simulation and dynamical 
theory. In practice, however, we do not have adequate in- 
formation about strengths, growth rates or duty cycles of 
spiral features. So empirical parametrization has so far been 



the only possible approach (Jenkins & Binney 1990; Jenkin; 



1992|)J 

Lacey and Ostriker (1985) proposed that the Galactic 
dark halo might consist of massive black holes (BHs) which 
would heat the disc as they pass through it. The total heat- 
ing produced by their proposed 2 x lO'' M© black holes is 
consistent with observations, but the velocity ellipsoid was 
found to be quite round, which is inconsistent with observa- 
tions. Furthermore, if the dark halos of dwarf spiral galaxies 
were domi nated by similar B Hs their discs would be easily 
destroyed (Friese & al. 1995). Finally, BHs passing through 
the disc might be expected to accrete and reveal themselves 
as high proper motion X-ray sources, and no such objects 
have been seen. It seems worth noting though that the cen- 
tral black hole in the Milky Way demonstrates that the ac- 
cretion rate and emergent flux from a black hole of about 
t he desired rn ass (2.6 x 10^ M0) is far from well understood 
(Genzel 1998); the central black hole in the Milky Way is 
remarkably quiet in X-rays. Another way of detecting these 
putative black holes is by statistically studying large num- 
bers of orbits of nearby stars; this may be possible with very 
accurate distances and kinematics which will become avail- 
able for stars within some kiloparsecs from ESA's GAIA 
satellite. 

In order to avoid p roblem s related to the black hole sce- 
nario Carr and Lacey (1987) proposed dark clusters of less 



massive objects. The latest generation of cosmological sim- 
ulations of galaxy formation do indeed show quite clumpy 



dark ha los around Milky Way type spirals (Moore et al 
199S()~]rhe amount, size and distribution of these clumps 
is still very uncertain, mainly due to numerical resolution 
issues, but if they can survive in the inner regions of dark 
halos they may well be able to heat discs in a manner con- 
sistent with observations. 



With this in mind we have embarked on a new study 
of disc heating simulations, starting with examining heat- 
ing due to GMCs and BHs. The European Space Agency's 
Hipparcos satellite has provided a wealth of new data on 
the distances and kinematics of nearby stars, so that one 
can now construct very much improved measurements of 
ages and kinematics. The velocities of the stars are much 
improved through the excellent Hipparcos parallaxes, while 
the ages of stars can be much better estimated, particularly 
for stars near the disc main sequence turn-off, because of the 
greatly improved absolute magnitudes. 

Ideally, one would like to carry out a self-consistent N- 
body simulation of an entire disc galaxy, with and without 
the perturbing influences due to GMCs/BHs etc, but this is 
not yet possible even with the fastest computers. The com- 
putation would involve both the stars (at approx 1 M0) and 
the perturbers (for a range of masses, lO'^-lO^ Mo), imply- 
ing some 10^^ stars and thousands of GMCs to millions of 
BHs in the simulations. This is not yet feasible in an N-body 
simulation. In this paper we render the problem tractable by 
simulating the heating in a local patch of the Galactic disc 
consisting typically of few tens of thousands of (massless) 
tracer stars, embedded within a fixed Galactic background 
potential, through which massive perturbers move. By us- 
ing a local simulation method we can resolve very well the 
effects of perturbers on nearby stellar orbits, for compari- 
son with the locally obtained observati ons. A very similar 
method has been used by Fuchs et al. (1994) who studied 
disc heating due to GMCs and/or 5 x IO^Mq black holes. 
We are able to broadly reproduce the results of their code 
in this area, but we also explore a larger range of scenarios. 

From our simulations, we find that GMC heating is less 
effective than earlier thought and further that heating due 
to GMCs creates a flatter velocity ellipsoid than found in 
earlier studies. We find that BHs are not as effective in heat- 
ing the disc as earlier thought. However, very massive halo 
black holes {Mbh = 1 x 10^ Mq) can heat the disc up to 
the observed amount, but the results of the heating rate 
and the ratio of the velocity dispersions are only marginally 
consistent with observations. A well selected combination of 
BHs and GMCs can give the observed heating in every com- 
ponent of the velocity dispersion within the observational 
error limits, but while interesting this solution is probably 
ad hoc, as the simulations do not yet model the effects of 
other perturbation sources, such as spiral waves. This will 
be examined in future work. 

This paper is organised as follows. In section 2 we re- 
view pre- and post-Hipparcos work on the observational de- 
termination of the age-velocity relation (AVR) for the lo- 
cal disc. We find the AVR is much less well measured than 
has previously been assumed. Observational issues thus still 
fundamentally limit our ability to discriminate between disc 
heating mechanisms. In section 3 we describe briefly our 
simulation method and numerical accuracy. In sections 4, 5, 
and 6 we demonstrate how massive perturbes (GMCs and 
halo BHs) heat up the stellar disc. We draw our conclusions 
in section 7. 
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2 OBSERVATIONS 

The Hipparcos satellite measured parallaxes and space mo- 
tions for a complete sky coverage of about 60,000 stars, 
with a further 60,000 stars assembled from a variety of pre- 
existing sources. As a result, the kinematics of stars used to 
construct the relationship between mean stellar age and the 
components of the velocity dispersion (the Age- Velocity re- 
lation, AVR) have been greatly improved. We briefly review 
here the state of the observational AVR both pre-Hipparcos 
and post-Hipparcos. 



2.1 Observational Age- Velocity- Relation 

Figure 1(a), shows several determinations of the AVR 
from the literature. We divide the measurements into pre- 
Hipparcos and post-Hipparcos determinations. 

Each AVR presented and the simulation results has 
been fit with the following form 



T 



(1) 



where a{t) is the total velocity dispersion as a function 
of time, t, (Jo is the initial velocity dispersion at f = 0, r 
i s a constant with unit of tim e, and a is the heating index 
(Wielen 1977; Villumsen 1985). Our motivation here is to fit 
the data with a simple analytical law which can also be fit 
to the simulations. 



2.1.1 Pre-Hipparcos determinations 

Firstly, consider the pre-Hipparcos measurements of the 
AVR, whic h are shown in panel ai of Figure 1. Edvards- 
son et al's (1993) data set of accurately determine d age s for 
189 F and G stars is shown by crosses. Wielen's (1977) de- 
termination, based mainly on K and M dwarfs which have 
been ranked by age based on chromospheric emission line 
measurements, is shown by square s. An improved version 
of the data set ( Fuchs fc al. 2000[ ) is shown by triangles. 



These latter data sets are mostly M dwarf types, for which 
the kinematical and parallax measurements have only been 
marginally improved by Hipparcos, and they remain robust 
in the post-Hipparcos epoch. 

We have fit each data set to Eqn |^, using a non-linear 
least squares gradient-expansion algorithm (using the IDL 
software package). All data points were given equal weight. 
For these three data sets, we find heating indices a in the 
range 0.47 to 0.61 (see Table 1). 

Our motivation is to fit the same form of the heating 
law to all the available data and to all the simulations. We 
should note that the heating indices a we derive are not quite 
the same as those derived by the authors above, who apply 
additional physical constraints to the fitting (e.g. the nature 
of the stellar orbital diffusion, and the weight to be applied 
to the initial value of the velocity dispersion). Nevertheless, 
the fitted indices are similar. 



2.1.2 Post-Hipparcos determinations 

We now consider post-Hipparcos AVR deter minat ions (Fig- 
ure 1(02)). Circles show data from Holmberg (2001), for 1486 



Table 1. Values of the heating index obtained by fitting data 
in the literature. The first five rows show our fits, and the last 
two show the author's own reported values. Note that Holmberg 
(2001) removes the effect of higher velocity thick disc stars when 
making his (more reliable) estimate of the heating index, whereas 
we leave them in (in order to estimate the bias they impose). 
There is a very considerable range in the derived values of a, 
which far exceeds the reported (internal) error estimates; a is 
therefore very likely to be dominated by systematic errors in the 
methods used to estimate stellar ages. 

Study a 

R. Wielen (1977) 0.61 ± 0.04 

Fuchs et al. (2000) 0.59 ± 0.04 

Edvardsson et al.(1993) 0.47 ± 0.002 

H. Rocha-Pinto 0.26 ± 0.008 

Holmberg (2001) All stars 0.45 ± 0.04 
Holmberg (2001) Thick Disc removed 0.33 ± 0.03 

Binney et al. (2000) 0.33 ± 0.03 



Hipparcos F and G type stars for which ages can be deter- 
mined with an accuracy of about 30% (of the age), by fit- 
ting absolute magnitude and colours to stellar isochrones. 
We show by diamonds our analysis of data kindly provided 
in advance of publication by H. Rocha-Pinto, who has de- 
termined ages of 425 mainly F, G and K stars using a care- 
f uUy calibrated technique bas ed on chromospheric emission 
(Rocha-Pinto fc Maciel 199^). For his stars we have mea- 



sured the total velocity dispersion as a function of age for 
the objects with [Fe/H]> —0.5, being careful to exclude sev- 
eral high velocity outliers. 

When fitting the data from H. Rocha-Pinto we yield a 
heating index a = 0.26 with a very small (formal) error (see 
Table hi). However, the value of the heating index is sensitive 
to the fitting procedure. For example, if we had used bins 
of 60 stars instead of 40-star-bins, we would have obtained 
a = 0.29±0.02. Because most of the stars in the sample were 
young we binned the data so that there were equal number 
of stars in each bin rather than using bins of equal width in 
age. 

The data from Holmberg was fitted in similar man- 
ner: the objects with [Fe/H] < —0.5 were removed from 
the data set and equal number bin were used. We obtained 
a = 0.45 ± 0.04 for the heating index. This is close to Holm- 
berg's own results (a — 0.46±0.02) for all his stars. However, 
Holmberg extends his fitting procedure by separating the 
stellar data into thin and thick disc components by a maxi- 
mum likelihood method, deriving a best fitting AVR for the 
thin disc for which the h eating index is a = 0.33 ± 0.03. 

Binney et al. (2000) have used a different approach to 
constrain the star forming history and velocity dispersion 
evolution in the Solar neighbourhood. They have used pho- 
tometrically complete sample of 11,865 Hipparcos stars. By 
combining local colour-magnitude diagram with the Padua 
isochrones and, crucially, kinematic information, they have 
obtained a heating index of a = 0.33 ±0.03 and a disc age of 
11.1 ±0.8 Gyr. The value obtained is consistent with Holm- 
berg's, for a smaller set of stars drawn from the same source 
(Holmberg works with individually determined stellar ages, 
and limits himself only to those for which the age can be 
determined quite accurately). 
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Figure 1. The observed total stellar velocity dispersion and ratio of velocity dispersions for nearby disc stars, as obtained by different 
authors. Best fits to the heating laws as parameterised by Eqn. f are also shown. Panel (ai): the solid line is the fit to Fuchs et al. data 
(triangles), the dashed line is the fit to Wielen's data (squares), and the dot-dashed line is the fit to Edvarsson et al. data (crosses). 
Panel (a2): the solid line is the fit to Rocha-Pinto's data (diamonds) and the dashed line is the fit to Holmberg's data (circles). Panels 
(bl) and (62) show the observed ratios of the vertical to the radial velocity dispersion for the same data sets and using the same symbols. 



While the values of a obtained by Binney et al. (200C) 



and Holmberg (2001) are consistent, they are inconsistent 



with the value we derive from Rocho-Pinto's data set, which 
uses an independent method to estimate stellar age, and in- 
consistent with the heating index derived for the K and M 
dwar fs, which have been age ranked via chro mospheric emis- 
sion (iFuchs fc al. 2OOOI). Fuchs et al. (|200(]|) show that the 



discrepency between values for the heating index of 0.3 and 
0.5 is actually rather mild (see their Fig. 1), when plotted 
in the age-velocity dispersion plane. 

The data sets give significantly different values for the 
heating index, ranging from 0.3 to 0.5. The error in the 
ages does give rise to a small systematic error in a; we have 
performed Monte-Carlo simulations in which we change the 
ages of the observed stars by a Gaussian distributed random 
fractional error, and recomputed the fit for a. Random age 
errors of order 20% typically increase a slightly, typically 
by less than 0.05. A 20% error in stellar ages is realistic in 
the context of determining the heating index in the AVR. 
For example, Holmberg uses only stars with a fractional er- 
ror of typically 20-30% (and always less than 50%), while 
Rocha-Pinto (private communication) reports fractional er- 
rors of order 30%. However, for the youngest stars in the 
sample, the errors may be significantly larger, of order 0.5 
in log(age). Monte-Carlo simulations show that this is cer- 
tainly enough to reduce the heating index from 0.5 to 0.3. 
Furthermore, the determination of a is sensitive to system- 
atic errors in the ages of the younger stars. We have per- 
formed Monte-Carlo simulations which show that changes 
in the value of a of order 0.1 can be obtained if, for exam- 
ple, there is a systematic error of 50% in the assigned ages 
of stars younger than a few Gyr (as might occur as a result 
of systematic errors in isochrone colours). There are thus 



several mechanisms which could produce systematic errors 
in the heating index as derived from observations. 

We are reluctant to arbitrarily assign any of these data 
sets greater credence than any other; they are all very care- 
fully constructed and they use different techniques, some 
independent, to assign ages to the stars. Instead we move 
on to consider the ratio of the velocity dispersions, rather 
than the vertical velocity dispersion alone. This quantity 
turns out to be much less sensitive to the ages of the stars 
and turns out to be very useful in comparing the data to the 
simulations. 



2.2 Ratio of Velocity Dispersions 

In panels (61) and (62) of Figure 1 the ratio of the vertical to 
the radial velocity dispersion Gzjo^ is plotted for the same 
cases as in panels (ai) and (02) and using te same symbols. 
The scatter seems to be larger amongst the pre-Hipparcos 
observations (panel (&i)) than in the post-Hipparcos data 
(panel (&2)). The latter observations indicate that the ratio 
is Ozjox — 0.5 and the former observations corro borat e this. 
Using the Hipparcos data, Dehnen and Binney ( 1998| ) have 
determined the velocity dispersion evolution as a function 
of colour index B — V. Their results also indicate that the 
ratio az/o-x tends to ~ 0.5. 



2.3 Comparing observations and simulations 

The observations show that the heating index lies between 
a = 0.3 and 0.5: in other words a is not yet well constrained 
by the observations. This was a major uncertainty in com- 
paring the results of the simulations and observations. On 
the other hand, the ratio az/ax is quite well measured at 
0.5 ± 0.1. The ratio of vertical velocity dispersion to radial 
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velocity dispersion thus turned out to be the best discrimi- 
nator when theoretical models or numerical simulation are 
compared with observations (Sections 4, 5, and 6). 



3 SIMULATION METHOD 

The local simulation method was first us ed in planetary 
ring dynamics by Wisdom and Tremaine (1988) (see also 



Salo 1 995) and was first applied to disc galaxies by Toomre 
(1990). In the method all calculations are restricted to a 



small co-moving box within the Galactic disc. A typical com- 
puter run simulates a "patch" of the disc, being 1 x 1 or 2 x 2 
kpc square, and 2 to 4 kpc in vertical extent, containing of 
order 10'* tracer stars in the disc as well as any perturbers 
(GMCs or BHs). 

The coordinates of stars in the region are referred to a 
point which moves on a circular orbit around the Galactic 
center at a distance r. A rotating Cartesian coordinate sys- 
tem is placed at the reference position, the a;-axis pointing 
radially outward (away from the Galactic center), the y-axis 
in the direction of orbital motion, and the z-axis in the direc- 
tion normal to the Galactic plane. The orbits of stars in the 
patch are computed by integrating the linea rized equations 
of motion ( |Hill 187^ [Julian fc Toomre 196^ ): 



y + 2Q.X = Fy 



z + u z 



F, 



(2) 



where Q is the orbital frequency, and k and v are the 
epicyclic frequency and the vertical frequency of motion, 
r espectively. The numerical integrator is the RK4 routine 
(Press & al. 1992). We adopt Solar neighbourhood (r = 8 
25.9 kms~ ^/kpc, k = 36.0 kms~^ / kpc. 

The 



kpc) values of fl 

and v = 98.7 km s~^/ kpc (Binney 



Tremaine 1987 



tracer stars feel the Galactic potential and the forces due to 
perturbers (GMCs or black holes). In order to avoid numer- 
ical artefacts, the gravitational forces {Fx, Fy,F^) are calcu- 
lated so that Newton's third law is satisfied, i.e. a perturber 
also feels the gravitational force from tracer stars. 



3.1 Numerical accuracy 

All the calculations are performed in double precision. In 
an unperturbed test run with an integration step size At = 
T'or6/800 the relative error in the stellar vertical kinetic en- 
ergy was observed to be \/:\Ez/E^{Q)\ < 4 x 10"^ at the 
end of the simulation run {I-end ~ 50To,.b ^ 12 Gyr). Here 
Torb — 237.23 X 10^ years, is the orbital time around the 
Galactic center. During the run the vertical component of 
the system angular momentum relative to the reference po- 
sition should remain zero, and is conserved to better than 



lA/zl/y /|(0) < 3 X 10 in which we have compared the 
change in the system angular momentum to the mean initial 
angular momentum of the stars. 

In the production runs we have used a still shorter inte- 
gration step size: from At = Torb/lOOO to At = Torb/ieOOO 
depending on the size and mass of the perturber. We have 
used such a short time step because the numerical accuracy 
of the orbits is dominated by the close encounters with the 



massive perturbers. Most of the numerical error is due to 
these instances. The time step would be even shorter if we 
were interested in the orbits, but it is actually the statisti- 
cal heating of many orbits which concerns us, rather than 
maximum accuracy for particular orbits. The time step was 
determined by simply testing how a shorter and shorter time 
step affects the evolution of the velocity dispersion. Because 
the integration error creates artificial noise in the system, 
we just need to search for the largest step size for which the 
velocity dispersion of the tracer stars does not increase due 
to step size in a simulation with massive perturbers. The 
optimal step size will naturally be different for a point mass 
like black holes and for an extended object like a GMC. 

An important potential source for systematic error is 
in the calculation of gravitational force. In practice, because 
we are using a local simulation box, the calculation of the 
gravitational interaction has to be cut off at some distance. 
Furthermore, in a local simulation method this limit should 
be at most half of the box size, in order to avoid comput- 
ing the gravitational force from a perturber more than once. 
We need to choose a box size for which the heating of the 
orbits by distant encounters outside the cut-off-radius is in- 
significant compared to closer encounters within the cut-off- 
radius. We have performed sets of simulations in order to 
test for the appropriate patch size. Some examples will il- 
lustrate this: for a 1 x 1 kpc^ patch and black holes with 
Mbh = 1 X 10® Mq we find the total stellar velocity dis- 
persion after 50 Solar orbital periods rises to atot = 22.7 
kms~^; changing to a patch of 2 x 2kpc'^, and running the 
simulation again with a larger cutoff length for gravity, we 
derive a total velocity dispersion of 22.1 kms"'^ at the end of 
the run. The difference between the runs is within the sim- 
ulation noise level and we conclude that the effects of long 
range forces are simulated with good accuracy. For GMCs 
of similar mass the patch size 1x1 kpc'^ was observed to be 
too small, and minimum size 2x2 kpc^ had to be used. 



3.2 Boundary conditions and the tracer stars 

Tracer stars move on epicycles in the patch, and can po- 
tentially cross the patch's boundaries. Periodic boundary 
conditions are used to recover these stars, which enter from 
the opposite boundary with a correction made for Galac- 
t ic sh ear, as described in detail by Wisdom and Tremaine 
The box in which the calculations are performed is 



( [1988D . 



thus effectively surrounded by virtual boxes which shear due 
to Galactic rotation, as illustrated in Figure 2. We typically 
use boxes which are 1x1 kpc or 2 x 2 kpc in the Galactic 
plane, and extend 2 or 4 kpc above and below it. 

The linearized equations are valid as long as and 
\z\ <C r. For all our simulations the velocity dispersions re- 
main low and this condition is met. 



3.3 Massive Perturbers 

We simulate two types of massive perturber: Giant Molecu- 
lar Clouds (GMCs) and Black Holes (BHs). The GMCs are 
placed in the Galactic plane on near-circular orbits, which 
are then integrated in the same manner as the tracer stars 
(but without the effect of the other perturbers, as this would 
unnecessarily slow the simulation and heat up the GMC 
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Figure 2. The simulation box (solid line) and its surrounding vir- 
tual boxes (dashed lines) are illustrated. Gravitational forces on a 
given star (cross) are calculated from the gravitational perturbers 
whose nearest image lies within the distance Rg (denoted by the 
circle). If a star (asterisk) leaves the simulation cell, its image will 
enter the cell in a way defined by the boundary conditions. 



population). The same periodic boundary conditions which 
replaces those stars which leave the box are also applied to 
the GMCs. 

The BiHs are on orbits characteristic of the dark halo, 
with high velocities relative to the disc stars. BiHs move 
through the box on essentially straight line orbits. They 
obey the same periodic boundary conditions as the tracer 
stars in the radial (a;-) and azimuthal (y-) direction, but not 
in the vertical (2-) direction. The border in the z-direction 
is given special status, by setting it to be far enough from 
the disc midplane, in order that no forces need to be calcu- 
lated across it. In practice one must take care that the box 
is large enough in the vertical direction, so that the stars 
never get close to the vertical boundary. When a halo black 
hole goes through the vertical border, it is removed from 
the simulation and a new one is randomly created from the 
parent velocity and space distribution. The number density 
of dark halo BHs entering and leaving the box is maintained 
at a constant level. 

The orbits of the tracer stars are computed by solving 
Equation ^ in the presence of massive perturbers in the disc 
or moving through the disc from the Galactic halo, by direct 
summation. The forces are calculated on each tracer star 
for all perturbers which are within a radius Rg of the star 
(see Fig. 2). In practice we have always used the maximum 
distance allowed by the patch, i.e. the calculation radius is 
half of the box width. 



4 DISC HEATING BY GIANT MOLECULAR 
CLOUDS 

A number of studies of disc heating via GMCs have been 
made, either via numerical simulations or via more theoret- 
ical approaches, e.g. by numerically integrating the Fokker- 
Planck equation. Our aim here is to compare our simulation 
results with previous work, and to compare with the post- 
Hipparcos observations of the AVR. 

Observations of the gas and dust in GMCs reveal a 
complex, fragmented structure that is often described as 
self-similar or fractal. We model our simulation GMCs by 
spheric al mass distribu t ions of uniform density. In many 
works dVillumsen 1983| ; |Villumsen 1985| ) the GMCs have 



been simulated as having a softened point potential (e.g. 
the Plummer model). However, when taking into account 
the complex structure of GMCs we have preferred using an 
homogeneous density distribution; we have thus modelled 
our GMCs as homogeneous spheres. The disc heating is not 
very sensitive to these differences in the modelling. 

Because most of the mass of interstellar matter is in the 
high end of the molecular cloud mass function, it is sufficient 
to simulate only the most massive GMCs, having typical 
masses of Mgmg = 1 x 10^ Mp) a nd diameters of dcMC = 
100 pc (|ScoviUe fc Sanders l"986|). Because GMCs do not 



interact with each other, the velocity dispersion of the GMC 
population is frozen during the simulation run. The observed 
velocity dispersion of the Galactic GMC population is so 
low that there must be some kind of physical mechanism for 
cooling it. There is some evidence that infalling gas into the 
Galactic disc cools down the interstellar gas. Whatever the 
physical mechanism, we can emulate this situation simply 
by allowing no interactions between the GMCs. 

We start the simulations with both the tracer stars and 
the GMCs having a low (cold) velocity dispersion, simulat- 
ing that the stars have just been born from the gas clouds, 
and share their kinematics. The values used are taken from 



Dehnen & Binney (1998). The total velocity dispersion a is 



decomposed into three components in the x, y, z directions 
of (axjCTyjO'z), with ax ~ 14 kms~^ and — 6 kms~^ 



{ay follows from the a^ and fi/ft-ratio). The initial veloc- 
ity distribution is Gaussian. The initial vertical structure 
is Gaussian. Strictly speaking the vertical structure should 
obey a sech'^-distribution in order to be isothermal, but in 
practice the sech^-distribution is very close to a Gaussian 
distribution in numerical simulations with a limited number 
of particles. 

In the GMC runs the simulation box size 2x2x4 kpc'' 
was found to be large enough for the calculation of grav- 
itational perturbations. After testing it was found that an 
integration step size At — Tort/ 1000 is sufficiently small for 
these runs, where Tort — 237.23 x 10® years. At the Solar 
radius the radial distribution of molecular clouds is falling 
of rapidly with increasing radius. In the vicinity of the So- 
lar neighbourhood the average density of molecular clouds 
i s approximately Egmc — 5 Mq / pc^ = 5 x 10® M0 / kpc'^ 
( ^coville fc Sanders 1986| ). 

The evolution of the stellar velocity dispersion is plot- 
ted in Figure 3(a) for four difi'erent assumed GMC number 
densities in the patch: Eqmc =2, 4, 8, 16 /kpc^. As most of 
the mass is at the high-mass end of the spectrum, the case 
Sgmc =4/ kpc'^ is closest to the present value for the Solar 
neighbourhood. 

We have fit the simulations with an age- velocity relation 
(AVR) of the form given in Eqn ^. The second curve, from 
bottom to top, in Figure 3(a) shows the fit with the most ap- 
propriate values for the Solar neighbourhood. The velocity 
dispersion is fitted only for the values atot > 30kms~^, be- 
cause the stellar heating by GMCs is known to obey a differ- 
ent power la w when the s tellar scale height is less than cloud 
scale height (Lacey 1984). Our fits are thus valid when the 



stellar population has a larger scale height than the GMC 
population. The velocity dispersion in the simulations is ob- 
served to grow as a{t) oc t°-^^ (see Table 2). This is slightly 
less than found in other studies {a{t) oc t"'^^ using differ- 



ent methods; i.e. analytical treatment in Lacey (1984), and 
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Figure 3. Panel (a). The evolution of total velocity dispersion a due to GMCs. Symbols are the simulation results and the curves 
represent least squares fits of the Eq. |^ to the points. Panel (b): The ratio of the vertical to the radial velocity dispersion shown as a 
function of the number density of GMCs, Sgmc> the simulation box. In both panels, from bottom to top, Sqmc = 2,4,8, 16kpc~^. 



"scaled" disc simulations in Villumsen (1985)). The least- 
square-fit-values for Q are shown in Table ^ individually for 
the total atot, the radial ax and the vertical component az 
of the velocity dispersion. The vertical velocity dispersion is 
observed to grow more efficiently than the velocity disper- 
sion in the horizontal pl ane. This tenden cy is in accoradance 
with the earlier results (Villumsen 1985), but our values for 



the Q are somewhat smaller. 

The ratio Ozjox of the vertical velocity dispersion to the 
radial velocity dispersion is plotted in a Figure 3(b). The ra- 
tio is clearly a function of the number density of the GMCs. 
At the end of the simulation this ratio varies between about 
0.45 and 0.60; i.e. the vertical heating is relatively most ef- 
fective when the GMC number density is highest. However, 
all the values fall within the observational limits az/^x — 
0.5 ±0.1. Coincidentally the mean value cjzjox — 0.5 is pro- 
duced by the case Eqmc ~ 4/kpc^ which is approximately 
the present GMC (of mass Mgmc = 1 x 10^ Mq) number 
density in the Solar neighbourhood. 

The ratio of the azimuthal component to the radial com- 
ponent depends on Q and k and is ay/ax — 0.7 for the Solar 
neighbourhood. This is also what is observed in the simula- 
tions. 

Our results seem to be in poor accordance with Lacey's 

(1984) results; he predicts az/(Jx — 0.8 for the present values 



of the epicyclic constants. He has also obtained a act for 
all the velocity components when the stellar scale height is 
larger than the GMC scale height. 

L acey g eneralized the work of Spitzer and Scwarzschild 

(1951; I953[ ) (who calculated the evolution of the stellar ve- 
locity distribution function by numerically integrating the 
Fokker-Planck equation) and calculated the evolution of 
all three components of the velocity dispersion by solv- 
ing the first moments of the Fokker-Plack equation. He as- 
sumed that the velocity distribution always remains a triax- 
ial Gaussian. Most importantly, the stellar orbits are as- 



sumed to be predominantly perturbed by many distant, 
weak encounters so that the evolution of the distribution 
function could be described by the diffusion equation. 

Lacey uses first-order epicyclic theory for the stellar or- 
bits (as we do here) and assumes the molecular clouds to be 
long lived (as here). However, for simplicity he adopts circu- 
lar GMC orbits, ignores Galactic shear and assumes that the 
interaction time between the stars and the GMCs is short. 

A simple test is to vary the GMC size: a simulation run 
with Egmc ~ 4/ kpc^ and doMC = 20 pc (which is similar to 
Lacey's cloud size) yielded a = 0.23 for the total velocity dis- 
persion. This goes in the expected direction as smaller cloud 
size will increase the strength of the gravitational perturba- 
tion in the close encounters. We also find that the vertical 
component of the velocity dispersion is increased relative to 
the radial one; we observed CTz/ctx — 0.57 at the end of the 
simulation. This partly explains the difference in the results. 
The remainder we ascribe to the assumptions in the analyti- 
cal method (i.e. circular GMC orbits, no Galactic shear and 
short star-CMC interactions). 

Numerical simulations by Villumsen (1985; 1983| ) have 
yielded dzjox — 0.6. He obtained similar results for the 
growt h of the velocity dispersion in the plane as Lacey 
(1984). On the other hand, for the vertical velocity disper- 



sion Villumsen obtained Oz <x t , which is significantly 
larger than the value we obtained (q = 0.26). 

Unfortunately, comparison of our results with Villum- 
sen's work is not straightforward because he analyses what 
one might term "scaled-up" simulations. Specifically, in his 
model, the exponential Galactic disc, with scale length 
1.75 kpc, contains 4000 GMCs with total mass Mtot = 
4x 10^ M0 outside 3 kpc radius, the individual GMCs having 
mass of A/gmc = 1 x 10® M© . This model is then scaled so 
that it could be simulated with fewer and heavier GMCs: 
the scaling being based on theoretical expectations (e.g. 
Lacey 1984) that in the diffusion process the efficiency is 
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Figure 4. Panel (a): the vertical density profile of the stellar disc at the start (filled circles) and at the end (open circles) of the simulation 
for a GMC number density of Eqmc = 16/kpc^. Panel (b): the vertical component of the velocity dispersion as a function of vertical 
height at the beginning (filled circles) and in the end (open circles) of the corresponding simulation. 



proportional to the total number of GMCs (A'^gmc) and 
to the square of the mass {Mq-^^q) of a single GMC, i.e. 
Nqmc X A^GMC = constant. In the actual simulations 200 
or 400 GMCs with mass Mgmc = 3.2 x 10*^ Mq have been 
used. The GMC size is the same as ours (dcMC ~ 100 pc), 
but a Plummer model has been used for the GMC potential 
and the cut-off radius for the gravitational interaction has 
been set at 2kpc. The numerical scheme is rather different 
from ours, but the equations of motion are integrated with a 
fourth order Runge-Kutta scheme as in our simulations. Vil- 
lumsen utilised an integration step size At = 3 x lO" years, 
corresponding to At — Torb/74: at the Solar radius. In hght 
of these many differences it is not surprising how difficult it 
is to compare our simulations with Villumsen's. 

One noticeable di fferen ce is, though, in the integration 
step size. Villumsen (1985) reported that relative numeri- 
cal error in energy was less than 0.01% in an unperturbed 
simulation. However, the errors in the numerical orbit inte- 
gration will be worst in close encounters between stars and 
massive GMCs. Based on our experience with our numerical 
scheme, a considerably larger step size than the one adopted 
by us induces artificial heating in the system. 

The vertical density distributions of the tracer stars are 
plotted in Figure 4(a) both at the beginning and the end 
of the simulation run with Egmc = le/kpc^. The initial 
density distribution (filled circles) is exactly Gaussian with 
scale height of Zo = 61 pc. At the end of the simulation run 
the vertical structure of the disc (open circles) is still almost 
Gaussian with an obviously much larger scale height {zo ~ 
230 pc) . Compared to a Gaussian distribution the high end 
tail is missing and the central density maximum is slightly 
suppressed, i.e. the profile is more 'boxy', but still very close 
to a Gaussian. 

The vertical velocity dispersion as a function of height is 
examined in Figure 4(b) for the same simulation. The initial 
velocity dispersion is isothermal (filled circles). At the end 



Table 2. Fitting of the GMC heating. The least-square-fit of the 
exponent a of the Eq. ^ is tabulated for the simulation runs along 
with the number density of GMCs. 



E 

[l/kpc2] 


a for crtot 


a for (7x 


a for (Tz 


2 


0.207 ±0.014 


0.198 ±0.019 


0.263 ±0.042 


4 


0.215 ± 0.010 


0.208 ±0.015 


0.265 ±0.031 


8 


0.217 ±0.008 


0.212 ±0.012 


0.263 ±0.024 


16 


0.217 ±0.007 


0.214 ±0.012 


0.252 ±0.015 


avg 


0.21 


0.21 


0.26 



of the simulation run, however, the velocity dispersion is 
observed to decrease as distance from the horiz ontal plane 
increases. Also numerical simulations by Jenkins (1992) have 



indicated that in the coeval disc the molecular cloud heating 
mechanism does not retain an isothermal velocity dispersion, 
but the velocity dispersion drops in higher altitudes as in our 
results. His simulations also yield similar results concerning 
the vertical density distribution; it is more 'boxy' than in 
the isothermal disc. 

To summarise the GMC section: our results rule out 
GMCs as the primary source of disc heating, because the 
GMC number density required to heat the disc by the 
correct amount is many times the observed GMC density 
in the Solar neighbourhood (and our results indicate that 
GMC heating is even less effective than earlier thought, i.e. 
a oc t^'^^). However, contrary to the earlier results, our sim- 
ulations indicate that the velocity ellipsoid is quite flat and 
consistent with the observational constraint on cTz/ax. 



Simulations of the heating of the Galactic stellar disc 9 



5 HEATING DUE TO HALO BLACK HOLES 

Another possibility for the heating of the Galactic disc would 
be due to gravitational interaction with massive black holes 
in the Galactic dark halo, w hich would make up pa rt or all 
of the Galaxy's dark matter (Lacey & Ostriker 1985). These 
authors estimate that black holes of mass Mbh « 2x10*^ M0 
orbiting inside the Galactic halo could explain the observed 
disc heating if they comprised the entire dark halo. 

One should note that there are clear difficulties with the 
black hole scenario. Based on estimates of the accretion rate 
of gas onto black holes which are passing through the disc, 
(McDowell 1985) rule out more massive black holes than 
Mbh ~ 1000 M0. Accretion should make such black holes 
directly observable as optical or infrared obje cts, bu t no such 
objects have been found. Lacey and Ostriker (1985) estimate 



that the accretion of the ISM onto the Mbh ~ 10^ 
halo black holes would create a number of X-ray sources 
which should have been observed if they existed. On the 
other hand, recent work on the Gal actic center b l ack hole, 
which has a mass of 2.6 x 10^ MQ(Genzel 1998; Narayan 



&Z. al. 1998), shows it to be a remarkably dim X-ray source 
despite the high density of gas in which it sits, which indi- 
cates that estimates of the X-ray flux of much nearer black 
holes of similar mass passing through the disc ISM may be 
overestimated by several magnitudes. 

Furthermore, this scenario has also purely dynamical 
difficulties. Very massive black holes in the Galactic halo 
might also disrupt globular clu sters in close encounters 
within the Hubble time. Moore (1993) estimates that the 
existence of the present population of low-mass globular 
clusters sets an upper limit of Mbh ~ 10' ' — 10 "* Mq for 



the black hole masses. Klessen and Burkert (1996) estimate 
the upper limit to be Mbh = 5 x lO** Mq . However, there 
are also some difficulties when deducing these upper limits. 
First of all, we do not know the initial mass function of glob- 
ular clusters. Secondly, as the evolution of a globular cluster 
is mainly affected by BH encounters within its small core 
region and as these events are rare, only few events during 
a Hubble time determine its fate. There will be thus large 
statistical variations in globular cluster lifetimes. 

Another limit on this scenario is introduce d by t he ef- 
fects of dynamical friction. Lacey and Ostriker (1985) have 
estimated that black holes of mass Mbh = 2 x 10*^ M© 
with initial orbital radii rnn < 2kpc would spiral into the 
Galactic center within t — 15 x 10® years. They assumed 
circular orbits for the black holes. The actual value for the 
spiralling radius depends on the background density distri- 
bution which is still not accurately known. 

Given these uncertainties, we attempt in this paper 
to establish constraints for the halo black hole population 
based purely on the kinematic evidence of heating in the 
Galactic disc, and ignore these other physical constraints. 

We test here the black hole scenario by simulating a 
patch of the local disc embedded in a dark Galactic halo 
composed all or partly of massive black holes. The black 
holes have been treated as point masses. We have simu- 
lated two different mass objects: Mbh = 1 x 10^ Mq and 
Mbh = 1 x 10^ Mq. We have choosen these masses be- 
cause smaller masses have an insignificant effect, and bigger 
masses are difficult to simulate in the present method {i.e. 
using a patch — the patch should be made so large that one 



would need to begin to take into account global properties 
of the Galaxy. In any case, black holes larger than 10* Mq 
can be ruled out because they would produce far too much 
heating, destroying the disc and/or destroying halo globular 
clusters ( |Carr 1994[ ). 

A softening of e = 1 pc has been used for the black holes 
in order to avoid numerical problems. 

The black holes are assumed to have zero net rotation 
around the Galactic center and a 1-D Gaussian velocity dis- 
persion of ax = cTy = cTz = 135kms~^. Black holes pass 
through the simulation box with randomly chosen velocities 
drawn from these kinematics. The stellar population has the 
same initial velocity dispersion as in the GMC simulation 
runs. 



5.1 1O''M0 Black Holes 

In the runs with Mbh = 1 x 10^ Mq, the patch size was set 
at 1 X 1 kpc^ in the horizontal plane. Different patch sizes 
were tested, but this size patch was found to be adequate. 
We have used 2 kpc for the vertical extent of the simulation 
box. The integration step size has to be small enough in 
order not to produce artificial heating. After some testing 
we found that At = ror6/8000 is small enough. 

We used 5000 stars in the patch and we investigated 
the effect of the number density of black holes on the disc 
heating. Figure 5(a) shows the evolution of the total velocity 
dispersion in five different cases: pbh ~ 4, 8, 16, 32, 64/ kpc''. 
These number densities correspond to mass densities 
PBH = 0.004, 0.008, 0.016,0.032, 0.064 MQ/pc^ respec- 
tively. For comparison, the mass density of the dar k halo 
in the Solar neighb ourhood is pn ~ O.OlMg/pc^ (Holm- 



berg fc Flynn 2000 ) , so the two higher number density runs 
are clearly unrealistic, and their purpose is mainly to illus- 
trate the effectiveness or non-effectiveness of this heating 
mechanism. 

It is already clear from Figure 5(a) that in order to 
account for the amount of disc heating that is observed in 
the Solar neighbourhood there should be an unrealistically 
large number of lO'' Mq black holes — the total mass of the 
black hole population would be larger than the mass of the 
whole Galactic halo. 

We have again used Eq. |l| in the least-square-fitting of 
the velocity dispersion. The resulting exponent a for the 
total, radial, and vertical velocity dispersions can be seen 
in Table 3. Our simulations verify the results of Lacey & 



Ostriker ( 1985 ) that the total velocity dispersion is propor- 
tional to (Tocr''^ for the heating due to the halo black holes. 

They also found that the vertical-radial-axis-ratio is 
az/^x ~ 0.67 in the non-selfgravitating case and a^/ax = 
0.55 in the self-gravitating case. For 10® Mq black holes we 
find that Ozjax ranges from 0.47 to 0.67, the highest ratio 
being for the highest black hole density (see Fig. 5b). The 
lowest value for a z I a x-tskAO is also closest to the observed 
value (Jz/(7x = 0.5 ± 0.1 

The development of the vertical density profile in the 
simulation with Mbh = 1 x 10® Mq and p = 32/kpc^ is 
shown in Figure 6(a). The density profile remains very close 
to a Gaussian, and has a half width of Zo « 350 pc. The 
vertical component of the velocity dispersion is plotted in 
Figure 6(b) as a function of vertical distance from the hor- 
izontal plane. The velocities remain isothermal up to high 
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Figure 5. Panel (a): the evolution of velocity dispersion a due to halo black holes of mass M = 1 X 10^ Mq. The curves represent 
simulation runs with black hole number densities pBH = 4, 8, 16, 32, 64/ kpc'^ from bottom to top. Panel (b): the evolution of the ratio 
of the vertical velocity dispersion to the radial velocity dispersion for the same simulations. 




Figure 6. Panel (a): The vertical density profile of the stellar disc is plotted at the beginning and end of the simulation. Panel (b): 
the vertical component of the velocity dispersion of the stellar disc is plotted as a function of height at the beginning and end of the 
simulation. The adopted black hole mass is Mbh = 1 x 10® Mq and the number density is pbh = 32/kpc^. 



altitudes, in this case up to about 700 pc at the end of the 
simulation. 



5.2 IO'^Mq Black Holes 

For the Mbh = 1 x 10^ M0 simulation runs a larger patch 
size was necessary. The used patch size was 2x2 kpc^ in the 
horizontal plane and 4kpc for the vertical size. The main 
reason for increasing the patch size is that we have to cal- 



culate gravity from larger distances with the more massive 
perturbers. 

The cutoff radius for calculating gravity was thus Rg = 
1 kpc. Testing showed that a smaller integration step size 
was also needed, and we adopted At = Tor(,/16000. 

For these simulations, 10, 000 stars were placed in the 
patch. We used black hole number densities of pbh = 
0.25, 0.5, 1.0/ kpc^, which correspond to mass densities 
PBH = 0.0025,0.005, 0.01 Mo/pc^, respectively. The evo- 
lution of the total velocity dispersion is plotted in Figure 
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Figure 7. Panel (a): the evolution of velocity dispersion a due to halo black holes of mass M = 1 X 10^ Mq. Panel (b): the evolution of 
the ratios of the vertical velocity dispersion to the radial velocity dispersion is plotted. The curves represent simulation runs with black 
hole number densities Pbh = 0.25,0.5, 1.0/ kpc^ from bottom to top. 



Table 3. Fitting of the disc heating index a due to the black 
holes of mass M = 1 X 10® Mq. 



Table 4. Fitting of the disc heating due to the black holes of 
mass Af = 1 X 10^ Mq. Also the standard deviation for the fitted 
parameter is given as an error estimate. 



p 

[l/kpc3] 


a for crtot 


a for (7x 


a for (Jz 


P 

[l/kpc3] 


a for atot 


a for ax 


a for (7z 










0.25 


0.52 ±0.04 


0.53 ±0.06 


0.50 ±0.04 


4 


0.529 ±0.014 


0.507 ±0.018 


0.597 ±0.030 


0.5 


0.48 ± 0.02 


0.48 ±0.02 


0.50 ±0.04 


8 


0.514 ±0.026 


0.524 ± 0.040 


0.498 ± 0.038 


1.0 


0.503 ± 0.008 


0.502 ±0.012 


0.509 ±0.016 


16 


0.535 ±0.057 


0.522 ± 0.082 


0.573 ±0.098 


avg 


0.50 


0.50 


0.50 


32 


0.53 ±0.13 


0.52 ±0.19 


0.56 ± 0.18 










64 


0.42 ±0.29 


0.42 ±0.40 


0.44 ±0.51 










avg 


0.51 


0.50 


0.53 


of the Galactic disc. However, if more weight is put on 



7(a). It is interesting to note that it would take about 8 bil- 
lion years to heat up the stellar disc from a ~ 18kms^^ 
to a ~ 80kms~^ if the total mass of the halo consisted of 
Mbh = 1 X 10^ M0 mass black holes. There are some am- 
biguity in the observations: some results indicate that the 
oldest stellar population would have a total velocity disper- 
sion only of (7 ~ 50 kms~^. In this case the required amount 
of heating would be produced if only half the halo mass were 
in the form of these black holes. The RMS-fits of a for differ- 
ent components of the velocity dispersion are found in Table 
4. 

In addition to the fact that these black holes broadly 
produce the right amount of disc heating, they also dis- 
tribute the heating in the three components consistently 
with observations. The (Tz/o"a:-ratios, plotted in a Figure 
7(b), range from 0.5 to 0.6, the highest black hole number 
density producing highest cTz/o-jj-ratio. Because of the am- 
biguities in the observed AVR, the massive halo black holes 
can not be ruled out as being responsible for the heating 



post-Hipparcos AVR determinations, then the value a — 0.5 
could be considered being too high. 

Our results are in mild disagreement with t he theoret- 



ical expectations of Lacey and Ostriker ( 1985 ). They es- 
timated that in order to heat the stellar disc from a « 
lOkms"^ up to (J ~ 80kms~^ in 15 Gyr the halo {p — 
0.01 Mq/ pc"^) should be comprised of black holes of mass 
Mbh — 2x 10*' Mq . Our simulations show that even Mbh = 
1 X 10^ M0 would not be enough to do the job. If the ini- 
tial velocity dispersion were around a ~ lOkms^^, the disc 
would be heated only to a ~ 60 kms^^ in 15 Gyr. (Note that 
Fig. 7(a) may be deceiving because a higher initial velocity 
dispersion, a ~ 18kms~^, was adopted). 



6 COMBINED EFFECT OF GMCS AND HALO 
OBJECTS 

In order to test the combined effects of the GMCs and halo 
black holes on the disc heating we performed two further sets 
of simulations. We know approximately the number density 
of the giant molecular clouds in the Solar neighbourhood, 
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Figure 8. Panel (a): the total velocity dispersion of the stellar disc. The combined effect of the GMC and black hole perturbations are 
seen with two different masses of the halo black holes: Mbh = 1 x 10^ Mq (filled circles) and Mbh = 1 x 10^ (open circles). Panel 
(b); the evolution of the ratios of the vertical velocity dispersion to the radial velocity dispersion. 



but none of the parameters related to the possible black hole 
population are directly known, except for the local density 
of the dark halo. 

The simulations seen in Figure 8 differ only by the num- 
ber density and mass of the halo black holes. In the run of 
Mbh = 1 X 10^ M0 the number density of the black holes 
is pBH = 8/kpc^ corresponding to pBH = 0.008 Mo/pc^. 
In this scenario most of the halo would thus consist of black 
holes. The combined effect of these black holes and GMCs 
would be enough to heat the disc to a = 50kms~^ in 12 
Gyr (filled circles in Figure 8(a)), which is somewhat less 
than wh at is a ctually observed (however, see also Dehnen & 
Binney (1998) who report a maximum velcocity dispersion 
a ~ 50 kms~^ for the oldest stellar population) . The growth 
of the velocity dispersion is now intermediate between the 
pure black hole and the GMC cases alone, a = 0.31. 

The ratio of the vertical and radial components of the 
velocity dispersion is plotted in Figure 8(b). In this case 
(filled circles) the ratio is around cr^/cTx — 0.53 after t ~ 10 
Gyr. This value fits well within observational error. 

The open circles in Figure 8 show the simulation with 
halo black holes of mass Mbh = 1 x 10^ M©, and a 
number density of pbh = 0.5/ kpc'^ which corresponds to 
Pbh ~ 0.005 M0/ pc'^. We can reproduce the disc heating 
very nicely with this setup, with the velocity dispersion ris- 
ing to a ~ 75 kms~^ over 12 Gyr, very similar to the obser- 
vations. The fitted heating law is within the observational 
limits at a = 0.42. The ratio of (Jz/(Jx is about 0.63, which 
is also within the observational error limits. 

It is interesting to note the high (Jz/gx ratios that the 
combined perturbations of BHs and GMCs produce. Why is 
the a^/a^ ratio not somewhere between the values obtained 
from the simulations in which BHs and GMCs act alone? 
The answer lies in the fact that (Jz/ctx ratio has a tendancy 
to grow when the number of perturbers is increased, as seen 
in both the GMC and BH simulation sets. When we increase 



the total number of perturbers, even if some of them are BHs 
and some GMCs, the az/cTx ratio is bound to increase. 



7 CONCLUSIONS 

There is firm observational evidence for the heating of the 
stellar disc with time. We have assembled the available post- 
Hipparcos samples of ages and kinematics for nearby stars 
and fit a heating law of the form a{t) = ao{l + t/T)'^ to 
the data, where a is the stellar velocity dispersion. The ob- 
servational data limit a to values between about 0.3 and 
0.6, i.e. despite the recent advances in measuring the kine- 
matics and ages of nearby stars due to the availablity of 
Hipparcos data, the age velocity relation in the nearby disc 
remains poorly determined. Fortunately, the ratio of the ver- 
tical component of the velocity dispersion to the radial com- 
ponent (Tzjox — 0.5 ± 0.1 is much better constrained by 
available data and can be used to test our simulations. 

We simulate disc heating by modelling a patch of the 
Galactic disc populated with tracer stars. We examine the 
effect of Giant Molecular Clouds (GMCs) travelling on disc 
orbits through the patch, and of massive black holes trav- 
elling through on plunging orbits characteristic of the dark 
halo. 

In our simulations of GMC heating we consistently find 
slightly lower values for a than earlier studies have found, 
which are based on either theoretical estimates of the dif- 
fusion rate of stellar orbits, or scaled simulations of the 
whole Galactic disc (with a small number of GMCs). We 
find a = 0.21 ± 0.02 for the evolution of the total velocity 
dispersion and a — 0.26 ±0.04 for the evolution of the verti- 
cal component of the velocity dispersion. The ratio a-z/ux of 
heating in the vertical and horizontal directions in the disc 
is found to depend on the adopted GMC density. 

The observed heating can not be explained by GMCs, 
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because even a four-fold GMC number density compared 
to the pres ent day value in the Sola r neighbourhood, E 
5Mq/pc'^ (Bcoville & Sanders 1986), does not heat up the 
disc enough. However, coincidentally we find that for the 
present day GMC number density the ratio Ozjox — 0.5 
which is consistent within the observational limits for stars 
in the local disc. 

We have also run simulations of the local disc heating 
being caused by massive (10^-10'^ Me) black holes of the 



dark h;io penetrating the disc and perturbing the stellar 
orbits. Uur simulations verity analytically obtained results 
in the literature, that the heating index is of order a = 0.5. 
Even if the entire halo consisted of Mbh — Ix 10^ Mq black 
holes, the resulting heating of the stellar disc would be much 
less than what is observed. However, if the whole halo were 
comprised of black holes of order Mbh = 1 x 10^ Mq the 
disc would be heated sufficiently. They would heat the disc 
in such a way that the ratio (Jz/ax would increase from 0.5 to 
0.6 which is still within the error limits of the observations. 

We have examined how the stellar disc would be heated 
by the combination of GMCs and halo black holes. Adopt- 
ing a GMC density consistent with the present observed 
value and a dark halo comprised of Mbh = 1 x 10^ Mq 
black holes would only heat the disc up to ct ~ 50kms~^ 
in 12 Gyr, which is inconsistent with most of the observa- 
tions. On the other hand, if half of the halo were made of 
Mbh = lx lO'' Mq black holes, they could with the GMCs 
heat the disc up to a ~ 75kms~^ in 12 Gyr, which is con- 
sistent with observations, but the heating would also push 
the ratio of the vertical velocity dispersion to the horizontal 
dispersion higher: a^jox — 0.63, which is still barely within 
the observational limits. In this case the heating index was 
observed to be a = 0.42. 

Taking into account the observational uncertainties in 
the stellar disc heating we can produce the right amount of 
the heating by a proper combination of massive halo black 
holes, Mbh = 1 x 10*^ Mq, and the right number density of 
GMCs in the Solar neighbourhood. Also the ratio a^jox and 
the heating index a fall within observational error limits. We 
do not regard this as a particularly satisfactory solution, 
because there are other potential causes of heating of stellar 
orbits, such as spiral arms, which we plan to investigate 
further in the simulation method presented here. 
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